Calculate metabolic index velocities for each species
Max Lindmark
06 March, 2023
Intro
Main analysis: This scripts reads predictions grids with predicted biomass densities (01_x) and metabolic index parameter estimates (02_x), calculates metabolic indexes by species and life stage and metabolic index velocities.
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(patchwork)
library(raster)
library(VoCC)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(viridis)
# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
theme_set(theme_plot())# Read predictions grids
# TODO
# This instead of ifelse...
# separate(id, c("species", "life_stage", "quarter", "var"), sep = "_", remove = FALSE) %>%
# mutate(species = str_to_title(species),
# life_stage = str_to_title(life_stage)) %>%
d <- read_csv("output/dat_pred.csv") %>%
separate(model, c("species", "life_stage"), sep = "_", remove = FALSE)
#> Rows: 7270592 Columns: 25
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): ices_rect, model
#> dbl (21): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div, yea...
#> lgl (2): sal_ct, sal_sc
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# mutate(species = NA,
# species = ifelse(model %in% c("cod_adult", "cod_juvenile"),
# "cod", species),
# species = ifelse(model %in% c("dab_adult", "dab_juvenile"),
# "dab", species),
# species = ifelse(model %in% c("flounder_adult", "flounder_juvenile"),
# "flounder", species),
# species = ifelse(model %in% c("plaice_adult", "plaice_juvenile"),
# "plaice", species))
# Read metabolic parameter estimates and left_join
mi_pars <- read_csv("output/mi_params.csv")
#> New names:
#> Rows: 4 Columns: 5
#> ── Column specification
#> ──────────────────────────────────────────────────────── Delimiter: "," chr
#> (1): species dbl (4): ...1, Ao_i, n, Eo
#> ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
#> Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> • `` -> `...1`
d <- left_join(d,
mi_pars %>% dplyr::select(species, Ao_i), by = "species")
#> left_join: added one column (Ao_i)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 7,270,592
#> > ===========
#> > rows total 7,270,592
# Read size csv to calculate the metabolic index
sizes <- read_csv("data/clean/sizes.csv") %>%
mutate(model = paste(species, name, sep = "_")) %>%
dplyr::select(model, B)
#> Rows: 8 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): species, name
#> dbl (2): mat_w, B
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'model' (character) with 8 unique values and 0% NA
d <- left_join(d, sizes, by = "model")
#> left_join: added one column (B)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 7,270,592
#> > ===========
#> > rows total 7,270,592d %>%
group_by(species) %>%
summarise(lwr = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.025)),
upr = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.975)))
#> group_by: one grouping variable (species)
#> summarise: now 4 rows and 3 columns, ungrouped
#> # A tibble: 4 × 3
#> species lwr upr
#> <chr> <dbl> <dbl>
#> 1 cod 2.12 11.7
#> 2 dab 1.28 12.1
#> 3 flounder 1.45 11.4
#> 4 plaice 2.14 13.3# First we need to convert oxygen in ml/L to PP
# Some conversions
# https://bluerobotics.com/learn/pressure-depth-calculator/
# First calculate pressure in pascal: p_tot = (r*g*h) + p_atm, then to kgPa and finally atm
rho = 1030 # density of water
p_atm = 1.01325 # atmospheric pressure
g = 9.807 # gravity of earth
d$po2 <- DO.unit.convert(d$oxy*(1/0.7),
DO.units.in = "mg/L",
DO.units.out = "PP",
bar.press = ((rho*g*d$depth + p_atm)*10^-5)*0.986923,
bar.units.in = "atm",
bar.units.out = "kpa",
temp.C = d$temp,
salinity = 10, # TODO: d$sal,
salinity.units = "pp.thou")
# Calculate the metabolic index for a given mass in the prediction grid a
# Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
kb <- 0.000086173324
t_ref <- 15 # arbitrary reference temperature
Eo <- unique(mi_pars$Eo)
n <- unique(mi_pars$n)
d <- d %>%
mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
phi = Ao_i*(B^n)*po2 * exp((Eo/kb)*inv_temp)) # Essington Eq. 1
#> mutate: new variable 'inv_temp' (double) with 908,489 unique values and 0% NA
#> new variable 'phi' (double) with 7,270,592 unique values and 0% NACalculate weighted quantiles of oxygen, temperature and phi by species and life stage
wms <- d %>%
pivot_longer(c(phi, oxy, temp), names_to = "env_var", values_to = "value") %>%
group_by(year, quarter, model) %>%
mutate(lwr = quantile(est, probs = 0.01),
upr = quantile(est, probs = 0.99)) %>% # trim quantiles
ungroup() %>%
filter(est > lwr & est < upr) %>%
group_by(env_var, year, quarter, model) %>%
summarise("Environment" = median(value),
"1st" = hutils::weighted_quantile(v = value, w = exp(est), p = c(0.25)),
"Weighted median" = hutils::weighted_quantile(v = value, w = exp(est), p = c(0.5)),
"3rd" = hutils::weighted_quantile(v = value, w = exp(est), p = c(0.75))) %>%
pivot_longer(c(Environment, `Weighted median`), names_to = "var", values_to = "Median") %>%
mutate(`1st` = ifelse(var == "Environment", NA, `1st`),
`3rd` = ifelse(var == "Environment", NA, `3rd`)) %>%
ungroup() %>%
mutate(life_stage = ifelse(grepl("juvenile", model), "Juvenile", "Adult")) %>%
mutate(species = gsub("_.*", "", model),
species = str_to_title(species))
#> pivot_longer: reorganized (oxy, temp, phi) into (env_var, value) [was 7270592x32, now 21811776x31]
#> group_by: 3 grouping variables (year, quarter, model)
#> mutate (grouped): new variable 'lwr' (double) with 448 unique values and 0% NA
#> new variable 'upr' (double) with 448 unique values and 0% NA
#> ungroup: no grouping variables
#> filter: removed 438,144 rows (2%), 21,373,632 rows remaining
#> group_by: 4 grouping variables (env_var, year, quarter, model)
#> summarise: now 1,344 rows and 8 columns, 3 group variables remaining (env_var, year, quarter)
#> pivot_longer: reorganized (Environment, Weighted median) into (var, Median) [was 1344x8, now 2688x8]
#> mutate (grouped): changed 1,344 values (50%) of '1st' (1344 new NA)
#> changed 1,344 values (50%) of '3rd' (1344 new NA)
#> ungroup: no grouping variables
#> mutate: new variable 'life_stage' (character) with 2 unique values and 0% NA
#> mutate: new variable 'species' (character) with 4 unique values and 0% NA# Group by quantile instead and use colors or linetypes
# Q1
pal <- brewer.pal(n = 3, name = "Spectral")
wms %>%
filter(env_var == "phi" & quarter == 1) %>%
pivot_longer(c(`1st`, Median, `3rd`)) %>%
ggplot(aes(year, value, linetype = var, color = name)) +
geom_line() +
scale_linetype_manual(values = c(2, 1)) +
#scale_color_brewer(palette = "Spectral", name = "Quantile") +
#scale_color_manual(values = pal, name = "Quantile") +
scale_color_manual(values = c("tomato2", "steelblue1", "grey20")) +
facet_grid(life_stage ~ species) +
labs(y = "Metabolic index (\u03C6)", x = "Year", linetype = "", color = "Quantile") +
theme(aspect.ratio = 1.2,
legend.position = c(0.1, 0.89),
legend.spacing.y = unit(0, 'cm'),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7))
#> filter: removed 2,240 rows (83%), 448 rows remaining
#> pivot_longer: reorganized (1st, 3rd, Median) into (name, value) [was 448x10, now 1344x9]
#> Warning: Removed 448 row(s) containing missing values (geom_path).
ggsave("figures/wm_mi_q1.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Removed 448 row(s) containing missing values (geom_path).
# Q4
wms %>%
filter(env_var == "phi" & quarter == 4) %>%
pivot_longer(c(`1st`, Median, `3rd`)) %>%
ggplot(aes(year, value, linetype = var, color = name)) +
geom_line() +
scale_linetype_manual(values = c(2, 1)) +
scale_color_manual(values = c("tomato2", "steelblue1", "grey20")) +
facet_grid(life_stage ~ species) +
labs(y = "Metabolic index (\u03C6)", x = "Year", linetype = "", color = "Quantile") +
theme(aspect.ratio = 1.2,
legend.position = c(0.1, 0.89),
legend.spacing.y = unit(0, 'cm'),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7))
#> filter: removed 2,240 rows (83%), 448 rows remaining
#> pivot_longer: reorganized (1st, 3rd, Median) into (name, value) [was 448x10, now 1344x9]
#> Warning: Removed 448 row(s) containing missing values (geom_path).
ggsave("figures/supp/wm_mi_q4.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Removed 448 row(s) containing missing values (geom_path).Plot slopes of weighted average quantiles over time
# Are there any "significant" trends over time?
yr_slopes <- wms %>%
filter(var == "Weighted median") %>%
mutate(id = paste(model, quarter, env_var, sep = "_")) %>%
split(.$id) %>%
purrr::map(~lm(Median ~ year, data = .x)) %>%
purrr::map_df(broom::tidy, .id = 'id') %>%
filter(term == "year") %>%
#separate(id_lw, c("year", "species"), sep = "_") %>%
mutate(lwr = estimate - 1.96*std.error,
upr = estimate + 1.96*std.error) %>%
separate(id, c("species", "life_stage", "quarter", "var"), sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) %>%
mutate(var = ifelse(var == "oxy", "Oxygen", var),
var = ifelse(var == "temp", "Temperature", var),
var = ifelse(var == "phi", "Phi", var))
#> filter: removed 1,344 rows (50%), 1,344 rows remaining
#> mutate: new variable 'id' (character) with 48 unique values and 0% NA
#> filter: removed 48 rows (50%), 48 rows remaining
#> mutate: new variable 'lwr' (double) with 48 unique values and 0% NA
#> new variable 'upr' (double) with 48 unique values and 0% NA
#> mutate: changed 48 values (100%) of 'species' (0 new NA)
#> changed 48 values (100%) of 'life_stage' (0 new NA)
#> mutate: changed 48 values (100%) of 'var' (0 new NA)
# Plot year slopes!
ggplot(yr_slopes, aes(x = factor(var, level = c('Oxygen', 'Temperature', 'Phi')),
estimate, ymin = lwr, ymax = upr, color = species, shape = quarter)) +
geom_hline(yintercept = 0, linetype = 2, color = "grey60") +
geom_point(position = position_dodge(width = 0.6)) +
geom_errorbar(width = 0, position = position_dodge(width = 0.6)) +
scale_color_brewer(palette = "Dark2") +
facet_wrap(~life_stage) +
labs(x = "", y = "Year slope", shape = "Quarter", color = "Species") +
guides(color = guide_legend(order = 1)) +
scale_x_discrete(labels = c('Phi' = expression(phi))) +
theme(aspect.ratio = 1.5,
legend.position = c(0.07, 0.88),
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7))
ggsave("figures/yr_slopes_mi.pdf", width = 17, height = 17, units = "cm")
# wms %>%
# filter(env_var == "oxy" & var == "Environment" & quarter == 1) %>%
# ggplot(aes(year, Median, color = life_stage, linetype = species)) +
# geom_line()
#
# wms %>%
# filter(env_var == "temp" & var == "Environment" & quarter == 1) %>%
# ggplot(aes(year, Median, color = life_stage, linetype = species)) +
# geom_line()
#
# wms %>%
# filter(env_var == "phi" & var == "Environment" & quarter == 1) %>%
# ggplot(aes(year, Median, color = life_stage, linetype = species)) +
# geom_line()Calculate phi velocities!
# Loop through year, quarter, species and life stage, make a dataframe of phi-velocities
d2 <- d %>%
dplyr::select(year, quarter, phi, model, X, Y) %>%
mutate(id = paste(quarter, model, sep = "_"))
#> mutate: new variable 'id' (character) with 16 unique values and 0% NA
velo_list <- list()
for(i in unique(d2$id)) {
d_sub <- filter(d2, id == i)
# Convert pred grid to raster stack
# https://gist.github.com/statnmap/a3eb564d06dee47f2d83eceac4212881
d_sub_rs <- d_sub %>%
group_split(year) %>%
map(~rasterFromXYZ(.x[,c("X", "Y", "phi")])) %>%
stack()
#plot(d_sub_rs)
# Temporal trend
vt2 <- tempTrend(d_sub_rs, th = 10) # th = Integer minimum number of observations in the series needed to calculate the trend at each cell.
# Spatial gradient
vg2 <- spatGrad(d_sub_rs, projected = TRUE)
# Climate velocity
gv2 <- gVoCC(vt2, vg2)
# (so that we can use ggplot)
points <- rasterToPoints(gv2)
# Make the points a dataframe for ggplot
df <- data.frame(points) %>%
rename(X = x, Y = y) %>%
mutate(id = i) %>%
separate(id, c("quarter", "species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
# Save!
velo_list[[i]] <- df
}
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
phi_velo_df <- dplyr::bind_rows(velo_list)Calculate biotic velocities!
# Loop through year, quarter, species and life stage, make a dataframe of phi-velocities
d2 <- d %>%
dplyr::select(year, quarter, est, model, X, Y) %>%
mutate(id = paste(quarter, model, sep = "_"),
exp_est = exp(est))
#> mutate: new variable 'id' (character) with 16 unique values and 0% NA
#> new variable 'exp_est' (double) with 7,270,592 unique values and 0% NA
velo_list <- list()
for(i in unique(d2$id)) {
d_sub <- filter(d2, id == i)
# Convert pred grid to raster stack
# https://gist.github.com/statnmap/a3eb564d06dee47f2d83eceac4212881
d_sub_rs <- d_sub %>%
group_split(year) %>%
map(~rasterFromXYZ(.x[,c("X", "Y", "exp_est")])) %>%
stack()
#plot(d_sub_rs)
# Temporal trend
vt2 <- tempTrend(d_sub_rs, th = 10) # th = Integer minimum number of observations in the series needed to calculate the trend at each cell.
# Spatial gradient
vg2 <- spatGrad(d_sub_rs, projected = TRUE)
# Climate velocity
gv2 <- gVoCC(vt2, vg2)
# (so that we can use ggplot)
points <- rasterToPoints(gv2)
# Make the points a dataframe for ggplot
df <- data.frame(points) %>%
rename(X = x, Y = y) %>%
mutate(id = i) %>%
separate(id, c("quarter", "species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
# Save!
velo_list[[i]] <- df
}
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
#> filter: removed 6,816,180 rows (94%), 454,412 rows remaining
#> rename: renamed 2 variables (X, Y)
#> mutate: new variable 'id' (character) with one unique value and 0% NA
#> mutate: changed 16,229 values (100%) of 'species' (0 new NA)
#> changed 16,229 values (100%) of 'life_stage' (0 new NA)
bio_velo_df <- dplyr::bind_rows(velo_list)Merge bio and phi velocity data frames, then add values from the prediction data frame (mean phi, mean biomass density)
# Cumulative filter instead
# d_test <- data.frame(y = round(rnorm(10, 10, 2))) %>%
# arrange(desc(y)) %>%
# mutate(cumsum = cumsum(y),
# sum = sum(y),
# perc = cumsum / sum)
#
# d_test
#
# d_test %>% filter(perc < 0.95)
# We want to regress the biotic velocity against the phi velocity, and we want to account for that biomass varies across grid cells, and we also want to see of the phi velocity has different effects depending on the average phi. We don't want to use all grid cells; some species and life stages occupy quite few cells. I will therefore take the mean biomass density across years for each grid cell, quarter and model, sort them and filter the grid cells that cumultively make up 95% of total biomass.
d_grid_keep <- d %>%
group_by(X, Y, quarter, model) %>%
summarise(grid_mean = mean(exp(est))) %>%
ungroup() %>%
group_by(quarter, model) %>%
arrange(desc(grid_mean)) %>%
mutate(cumsum = cumsum(grid_mean),
sum = sum(grid_mean),
perc = cumsum / sum) %>%
ungroup() %>%
filter(perc < 0.95) %>%
mutate(id = paste(X, Y, quarter, model, sep = "_"))
#> group_by: 4 grouping variables (X, Y, quarter, model)
#> summarise: now 259,664 rows and 5 columns, 3 group variables remaining (X, Y, quarter)
#> ungroup: no grouping variables
#> group_by: 2 grouping variables (quarter, model)
#> mutate (grouped): new variable 'cumsum' (double) with 259,099 unique values and 0% NA
#> new variable 'sum' (double) with 16 unique values and 0% NA
#> new variable 'perc' (double) with 204,648 unique values and 0% NA
#> ungroup: no grouping variables
#> filter: removed 165,298 rows (64%), 94,366 rows remaining
#> mutate: new variable 'id' (character) with 94,366 unique values and 0% NA
# Filter the main dataset to include these grid cells only, then calcualte mean phi and biomass in theses cells (to be used as covariates)
pred_grid_sum <- d %>%
mutate(id = paste(X, Y, quarter, model, sep = "_")) %>%
filter(id %in% unique(d_grid_keep$id)) %>%
group_by(id) %>%
summarise(mean_log_biomass = mean(est),
mean_phi = mean(phi)) %>%
ungroup() %>%
separate(id, c("X", "Y", "quarter", "species", "life_stage"),
sep = "_", remove = FALSE)
#> mutate: new variable 'id' (character) with 259,664 unique values and 0% NA
#> filter: removed 4,628,344 rows (64%), 2,642,248 rows remaining
#> group_by: one grouping variable (id)
#> summarise: now 94,366 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
# Now we can scale variables by model and quarter
pred_grid_sum <- pred_grid_sum %>%
group_by(quarter, species, life_stage) %>%
mutate(mean_log_biomass_sc = mean_log_biomass - mean(mean_log_biomass),
mean_phi_ct = mean_phi - mean(mean_phi)) %>%
ungroup() %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage)) %>%
mutate(velo_id = paste(X, Y, quarter, species, life_stage, sep = "_")) %>%
dplyr::select(velo_id,
mean_log_biomass, mean_log_biomass_sc,
mean_phi, mean_phi_ct)
#> group_by: 3 grouping variables (quarter, species, life_stage)
#> mutate (grouped): new variable 'mean_log_biomass_sc' (double) with 70,492 unique values and 0% NA
#> new variable 'mean_phi_ct' (double) with 94,366 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: changed 94,366 values (100%) of 'species' (0 new NA)
#> changed 94,366 values (100%) of 'life_stage' (0 new NA)
#> mutate: new variable 'velo_id' (character) with 94,366 unique values and 0% NA
# Join the biotic and phi velocities...
velo_df <- bio_velo_df %>%
mutate(id2 = paste(id, X, Y, sep = "_")) %>%
rename(bio_velo = voccMag) %>%
left_join(phi_velo_df %>%
mutate(id2 = paste(id, X, Y, sep = "_")) %>%
dplyr::select(id2, voccMag) %>%
rename(phi_velo = voccMag)) %>%
dplyr::select(-id2, -voccAng) %>%
drop_na(phi_velo) %>% # remove na and non-values
filter(is.finite(phi_velo))
#> mutate: new variable 'id2' (character) with 259,664 unique values and 0% NA
#> rename: renamed one variable (bio_velo)
#> mutate: new variable 'id2' (character) with 259,664 unique values and 0% NA
#> rename: renamed one variable (phi_velo)
#> Joining, by = "id2"
#> left_join: added one column (phi_velo)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 259,664
#> > =========
#> > rows total 259,664
#> drop_na: no rows removed
#> filter: removed 48 rows (<1%), 259,616 rows remaining
# ... then join in the grid cell variables by model + grid cell id (velo_id)
velo_df <- velo_df %>%
mutate(velo_id = paste(X, Y, quarter, species, life_stage, sep = "_")) %>%
left_join(pred_grid_sum, by = "velo_id") %>%
dplyr::select(-velo_id) %>%
drop_na(mean_log_biomass) # because we cropped the pred_grid_sum, these will be NA in velo_df
#> mutate: new variable 'velo_id' (character) with 259,616 unique values and 0% NA
#> left_join: added 4 columns (mean_log_biomass, mean_log_biomass_sc, mean_phi, mean_phi_ct)
#> > rows only in x 165,264
#> > rows only in y ( 14)
#> > matched rows 94,352
#> > =========
#> > rows total 259,616
#> drop_na: removed 165,264 rows (64%), 94,352 rows remaining
# Scale velocities before fitting
velo_df <- velo_df %>%
group_by(id) %>% # id is quarter + model (species + life_stage)
mutate(mean_phi_velo = mean(phi_velo),
sd_phi_velo = sd(phi_velo)) %>%
ungroup() %>%
mutate(phi_velo_sc = (phi_velo - mean_phi_velo)/sd_phi_velo)
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'mean_phi_velo' (double) with 16 unique values and 0% NA
#> new variable 'sd_phi_velo' (double) with 16 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: new variable 'phi_velo_sc' (double) with 94,352 unique values and 0% NA
# Filter extreme velocity values?
hist(velo_df$phi_velo)
hist(velo_df$bio_velo)
velo_df <- velo_df %>%
group_by(id) %>%
mutate(bio_velo_lwr = quantile(bio_velo, 0.005),
bio_velo_upr = quantile(bio_velo, 0.995),
phi_velo_lwr = quantile(phi_velo, 0.005),
phi_velo_upr = quantile(phi_velo, 0.995)) %>%
ungroup() %>%
mutate(bio_velo = ifelse(bio_velo < bio_velo_lwr, bio_velo_lwr, bio_velo),
bio_velo = ifelse(bio_velo > bio_velo_upr, bio_velo_upr, bio_velo),
phi_velo = ifelse(phi_velo < phi_velo_lwr, phi_velo_lwr, phi_velo),
phi_velo = ifelse(phi_velo > phi_velo_upr, phi_velo_upr, phi_velo))
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'bio_velo_lwr' (double) with 15 unique values and 0% NA
#> new variable 'bio_velo_upr' (double) with 16 unique values and 0% NA
#> new variable 'phi_velo_lwr' (double) with 16 unique values and 0% NA
#> new variable 'phi_velo_upr' (double) with 16 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: changed 960 values (1%) of 'bio_velo' (0 new NA)
#> changed 960 values (1%) of 'phi_velo' (0 new NA)
# Plot again
velo_df %>%
ggplot(aes(phi_velo)) +
geom_histogram() +
facet_wrap(~id, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
velo_df %>%
ggplot(aes(bio_velo)) +
geom_histogram() +
facet_wrap(~id, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Fit sdmTMB models of biomass velocities as a function of metabolic index velocities
# Make a for loop here..
coefs <- list()
pred_velo1 <- list()
pred_velo2 <- list()
for(i in unique(velo_df$id)){
dd <- filter(velo_df, id == i)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 150)
#mesh <- make_mesh(dd, c("X", "Y"), cutoff = 5)
plot(mesh)
m <- sdmTMB(
bio_velo ~ phi_velo_sc * mean_phi_ct + mean_log_biomass_sc,
mesh = mesh,
data = dd,
family = gaussian(),
spatial = "on",
spatiotemporal = "off",
control = sdmTMBcontrol(newton_loops = 1)
)
#sanity(m)
m1 <- run_extra_optimization(m)
sanity(m1)
# Save coefficients??
coefs[[i]] <- tidy(m1, conf.int = TRUE) %>% mutate(model = i)
# TODO: Should this really be within species? Guess so
# Now predict the biotic velocity for high and low temperature with a sd change in the climate velocity
low_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.025))
high_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.975))
sd_clim_v <- sd(dd$phi_velo)
mean_log_biomass_sc <- mean(dd$mean_log_biomass_sc)
# Predict
nd <- data.frame(phi_velo_sc = c(0, 1, 0, 1),
mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = 2),
mean_log_biomass_sc = rep(mean_log_biomass_sc, 4),
X = 0,
Y = 0)
pred <- predict(m1, newdata = nd, se_fit = TRUE)
delta_bio_v_low_phi <- pred$est[2] - pred$est[1]
delta_bio_v_high_phi <- pred$est[4] - pred$est[3]
nd2 <- data.frame(bio_vel_delta = c(delta_bio_v_low_phi, delta_bio_v_high_phi),
mean_clim = c("low_phi", "high_phi"),
model = i)
pred_velo1[[i]] <- nd2
# Now make a prediction across a range of velocities
# Again I'm predicting for the range of the specific model
nd3 <- data.frame(expand.grid(phi_velo_sc = seq(min(dd$phi_velo_sc), max(dd$phi_velo_sc), length.out = 50),
mean_phi_ct = c(low_mean_clim, high_mean_clim))) %>%
mutate(mean_log_biomass_sc = mean_log_biomass_sc,
X = 0,
Y = 0,
phi_velo = (phi_velo_sc*sd(dd$phi_velo)) + mean(dd$phi_velo))
#ggplot(nd3, aes(phi_velo, phi_velo_sc)) + geom_point() + geom_abline()
pred3 <- predict(m1, newdata = nd3, se_fit = TRUE)
pred3 <- pred3 %>%
mutate(lwr = est - 1.96*est_se,
upr = est + 1.96*est_se,
model = i)
pred_velo2[[i]] <- pred3
}
#> filter: removed 84,921 rows (90%), 9,431 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 84,921 rows (90%), 9,431 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 85,727 rows (91%), 8,625 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 85,727 rows (91%), 8,625 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 93,392 rows (99%), 960 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 93,392 rows (99%), 960 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 93,645 rows (99%), 707 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 93,645 rows (99%), 707 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 81,880 rows (87%), 12,472 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 81,880 rows (87%), 12,472 rows remaining
#> Warning: did not converge in 10 iterations
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 87,369 rows (93%), 6,983 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 87,369 rows (93%), 6,983 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 90,309 rows (96%), 4,043 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 90,309 rows (96%), 4,043 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 90,397 rows (96%), 3,955 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 90,397 rows (96%), 3,955 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_sc' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'phi_velo' (double) with 50 unique values and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'model' (character) with one unique value and 0% NA
dat_coefs <- dplyr::bind_rows(coefs) %>%
mutate(sig = "N",
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig),
sig = ifelse(estimate > 0 & conf.low > 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA
dat_pred_velo1 <- dplyr::bind_rows(pred_velo1) %>%
#left_join(dat_coefs_sig, by = "model") %>% # TODO fix this, need to filter specific params
separate(model, c("quarter", "species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 32 values (100%) of 'species' (0 new NA)
#> changed 32 values (100%) of 'life_stage' (0 new NA)
dat_pred_velo2 <- dplyr::bind_rows(pred_velo2) %>%
separate(model, c("quarter", "species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 1,600 values (100%) of 'species' (0 new NA)
#> changed 1,600 values (100%) of 'life_stage' (0 new NA)
write_csv(dat_pred_velo1, "output/dat_pred_velo1.csv")
write_csv(dat_pred_velo2, "output/dat_pred_velo2.csv")# Before plotting conditional effects and such, plot coefficients. We expect phi_velocity_sc to have a POSITIVE effect: the bigger the value, the LESS decline in ph. We expect the interaction to be NEGATIVE: The effect phi velocity declines as the mean phi is bigger!
# https://biologyforfun.wordpress.com/2014/04/08/interpreting-interaction-coefficient-in-r-part1-lm/
dat_coefs_sub <- dat_coefs %>%
filter(term %in% c("phi_velo_sc", "phi_velo_sc:mean_phi_ct")) %>%
dplyr::select(model, estimate, term) %>%
pivot_wider(names_from = term, values_from = estimate) %>%
mutate(Expected = ifelse(phi_velo_sc > 0, "Yes", "No")) %>%
mutate(Expected = ifelse(`phi_velo_sc:mean_phi_ct` < 0, "Yes", Expected)) %>%
dplyr::select(model, Expected)
#> filter: removed 48 rows (60%), 32 rows remaining
#> pivot_wider: reorganized (estimate, term) into (phi_velo_sc, phi_velo_sc:mean_phi_ct) [was 32x3, now 16x3]
#> mutate: new variable 'Expected' (character) with 2 unique values and 0% NA
#> mutate: changed one value (6%) of 'Expected' (0 new NA)
#dat_coefs %>% filter(term %in% c("phi_velo_sc", "phi_velo_sc:mean_phi_ct")) %>% arrange(desc(abs(estimate)))
dat_coefs_plot <- dat_coefs %>%
filter(term %in% c("phi_velo_sc", "phi_velo_sc:mean_phi_ct")) %>%
left_join(dat_coefs_sub, by = "model") %>%
mutate(Term = ifelse(term == "phi_velo_sc", "\u03C6 velocity", "\u03C6 velocity × mean \u03C6")) %>%
separate(model, c("Quarter", "species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(life_stage = str_to_title(life_stage),
species = str_to_title(species)) %>%
mutate(model = paste(life_stage, species))
#> filter: removed 48 rows (60%), 32 rows remaining
#> left_join: added one column (Expected)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 32
#> > ====
#> > rows total 32
#> mutate: new variable 'Term' (character) with 2 unique values and 0% NA
#> mutate: changed 32 values (100%) of 'species' (0 new NA)
#> changed 32 values (100%) of 'life_stage' (0 new NA)
#> mutate: changed 32 values (100%) of 'model' (0 new NA)
ggplot(filter(dat_coefs_plot, sig == "Y"), aes(model, estimate, fill = Expected, color = Expected, shape = Term)) +
geom_hline(yintercept = 0, linetype = 2, alpha = 0.8) +
geom_point(size = 2, position = position_dodge(width = 0.4), color = "white", stroke = 0.4) +
geom_point(data = filter(dat_coefs_plot, sig == "N"),
size = 2, position = position_dodge(width = 0.4), color = "white", stroke = 0.4,
fill = "gray30") +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(width = 0.4), width = 0,
size = 0.6) +
geom_errorbar(data = filter(dat_coefs_plot, sig == "N"), size = 0.6,
aes(ymin = conf.low, ymax = conf.high), position = position_dodge(width = 0.4), width = 0,
color = "gray30") +
scale_shape_manual(values = c(23, 21)) +
scale_fill_manual(values = c("tomato2", "steelblue1")) +
scale_color_manual(values = c("tomato2", "steelblue1")) +
facet_wrap(~Quarter, ncol = 1) +
coord_flip() +
guides(shape = guide_legend(override.aes = list(shape = c(23, 21), color = "grey20", fill = "gray20"))) +
labs(y = "Estimate", x = "")
#> filter: removed 25 rows (78%), 7 rows remaining
#> filter: removed 7 rows (22%), 25 rows remaining
#> filter: removed 7 rows (22%), 25 rows remaining
ggsave("figures/mi_velo_pars.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
# Plot delta velocity
# Pivot back and forth to calculate difference between low and high phi (for sorting the plot)
# dat_pred_velo1 <- dat_pred_velo1 %>%
# pivot_wider(names_from = mean_clim, values_from = bio_vel_delta) %>%
# mutate(diff = low_phi - high_phi) %>%
# pivot_longer(c(low_phi, high_phi), names_to = "mean_clim", values_to = "bio_vel_delta")
#
# # Plot!
# dat_pred_velo1 %>%
# mutate(id2 = paste(species, life_stage)) %>%
# ggplot(aes(bio_vel_delta, fct_reorder(id2, diff), color = mean_clim, shape = sig)) +
# geom_point(size = 3, position = position_dodge(width = 0.1), alpha = 0.9) +
# facet_wrap(~quarter, ncol = 1) +
# scale_color_manual(values = c("steelblue", "tomato3")) +
# scale_shape_manual(values = c(16, 21)) +
# #scale_fill_manual(values = c("white", "grey10")) +
# geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
# NULL
dat_pred_velo2 %>%
mutate(id2 = paste(life_stage, species, paste0("Q", quarter))) %>%
filter(model %in% c("1_dab_juvenile")) %>%
mutate(mean_clim = ifelse(mean_phi_ct > 0, "High \u03C6", "Low \u03C6")) %>%
ggplot(aes(phi_velo, est, color = mean_clim, fill = mean_clim)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
geom_line(alpha = 0.8) +
labs(y = "Predicted Biotic velocity (km/year)",
x = "\u03C6 velocity (km/year)") +
scale_fill_manual(values = c("steelblue1", "tomato2"), name = "") +
scale_color_manual(values = c("steelblue1", "tomato2"), name = "") +
facet_wrap(~id2, scales = "free") +
theme(aspect.ratio = 1,
legend.position = c(0.92, 0.1),
legend.background = element_rect(fill = NA))
#> mutate: new variable 'id2' (character) with 16 unique values and 0% NA
#> filter: removed 1,500 rows (94%), 100 rows remaining
#> mutate: new variable 'mean_clim' (character) with 2 unique values and 0% NA
ggsave("figures/mi_velo_conditional.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)Extra plots
# Species biomass filtered by the minimum biomass I use
# 4*2
d_plot <- d %>%
mutate(id = paste(X, Y, quarter, model, sep = "_")) %>%
filter(id %in% unique(d_grid_keep$id)) %>%
filter(quarter == 1 & year == 2010) %>%
separate(model, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: new variable 'id' (character) with 259,664 unique values and 0% NA
#> filter: removed 4,628,344 rows (64%), 2,642,248 rows remaining
#> filter: removed 2,595,065 rows (98%), 47,183 rows remaining
#> mutate: changed 47,183 values (100%) of 'species' (0 new NA)
#> changed 47,183 values (100%) of 'life_stage' (0 new NA)
plot_map +
geom_raster(data = d_plot, aes(X*1000, Y*1000, fill = est)) +
facet_grid(life_stage ~ species) +
scale_fill_viridis(alpha = 0.95, name = "Density [log(kg/km)]")
#> Warning: Removed 1309 rows containing missing values (geom_raster).
# plot_map +
# geom_raster(data = d_plot, aes(X*1000, Y*1000, fill = exp(est))) +
# facet_grid(life_stage ~ species) +
# scale_fill_viridis(alpha = 0.95, name = "Biomass density (kg/km)", trans = "sqrt")
ggsave("figures/biomass_pred.pdf", width = 20, height = 20, units = "cm", device = cairo_pdf)
#> Warning: Removed 1309 rows containing missing values (geom_raster).
# Plot velocities of metabolic index
# Not needed to save for each i, the fact that the absolute values are not the same does not matter because phi in the environment over time is driven by trends in temperature and oxygen also!
p_df <- velo_df %>% filter(species == "Cod" & life_stage == "Adult" & quarter == 1)
#> filter: removed 84,921 rows (90%), 9,431 rows remaining
plot_map +
geom_raster(data = p_df,
aes(X*1000, Y*1000, fill = phi_velo)) +
scale_fill_gradient2(name = expression("\u03D5 velocity (km"~""*year^-1*")"),
na.value = scales::muted("red"),
limits = c(quantile(p_df$voccMag, 0.01),
max(is.finite(p_df$phi_velo)))) +
geom_sf(size = 0.3) +
scale_fill_gradient2(name = "\u03C6 velocity (km/year)") +
theme(plot.title = element_text(size = 9),
legend.position = c(0.18, 0.9),
legend.direction = "horizontal",
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(0.6, "line"),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7),
legend.background = element_rect(fill = NA))
#> Warning: Unknown or uninitialised column: `voccMag`.
#> Scale for 'fill' is already present. Adding another scale for 'fill', which
#> will replace the existing scale.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
ggsave("figures/mi_velocity.pdf",
width = 14, height = 14, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
p_df2 <- d_plot %>% filter(species == "Cod" & life_stage == "Adult" & quarter == 1)
#> filter: removed 37,752 rows (80%), 9,431 rows remaining
p1 <- plot_map +
geom_raster(data = p_df2,
aes(X*1000, Y*1000, fill = oxy)) +
geom_sf(size = 0.3) +
scale_fill_viridis(option = "D", name = "Oxygen [ml/L]") +
theme(plot.title = element_text(size = 9),
legend.position = c(0.17, 0.9),
legend.direction = "horizontal",
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(0.6, "line"),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
legend.background = element_rect(fill = NA))
p2 <- plot_map +
geom_raster(data = p_df2,
aes(X*1000, Y*1000, fill = temp)) +
geom_sf(size = 0.3) +
scale_fill_viridis(option = "A", name = "Temperature [°C]") +
theme(plot.title = element_text(size = 9),
legend.position = c(0.2, 0.9),
legend.direction = "horizontal",
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(0.6, "line"),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
legend.background = element_rect(fill = NA))
p3 <- plot_map +
geom_raster(data = p_df2,
aes(X*1000, Y*1000, fill = phi)) +
geom_sf(size = 0.3) +
scale_fill_viridis(option = "mako", name = "\u03C6") +
theme(plot.title = element_text(size = 9),
legend.position = c(0.14, 0.9),
legend.direction = "horizontal",
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(0.6, "line"),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
legend.background = element_rect(fill = NA))
p1 + p2 + p3
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
ggsave("figures/spatial_vars.pdf",
width = 30, height = 20, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.# From VoCC tutorial
# devtools::install_github("JorGarMol/VoCC", dependencies = TRUE, build_vignettes = FALSE)
#
# library(VoCC)
#
# # monthly to annual averages
# library(repmis)
#
# # monthly to annual averages
# r <- sumSeries(HadiSST, p = "1960-01/2009-12", yr0 = "1955-01-01", l = nlayers(HadiSST), fun = function(x) colMeans(x, na.rm = TRUE), freqin = "months", freqout = "years")
#
# class(r)
# head(r)
#
# # temporal trend
# vt <- tempTrend(r, th = 10)
#
# # spatial gradient
# vg <- spatGrad(r, th = 0.0001, projected = FALSE)
#
# # climate velocity
# gv <- gVoCC(vt, vg)
#
# dv <- raster(gv)
#
# library(rasterVis)
# my.at <- seq(-50, 50, by = 5)
# p1 <- rasterVis::levelplot(gv[[1]], par.settings = BuRdTheme, at=my.at, main = 'Gradient-based vocc', margin = FALSE)